% Code based on Gertler, M. and P. Karadi (AJE: Macro, 2015),
% which is based on Mertens-Ravn, (AER, 2013) external instrument SVAR code.
% For the original code of Gertler and Karadi (2015) see:
% https://www.openicpsr.org/openicpsr/project/114082/version/V1/view

clear all; close all; clc;

% Import and load time-series data
AA_Import_Data;
load '../00_Data/DATASET' DATASET_VAR DATASET_FACTORS;

rng('default')
nboot  = 2000;      % Number of Bootstrap samples
clevel = 90;        % Bootstrap percentile shown
VAR.fontsize = 14;  % Font size in figures

% Setting up the specification
mtmp                            = 9
smpl_min_VAR_vec                = [2006 2];              % year month    
smpl_max_VAR_vec                = [2010 12];  			 % Choose sample
monpol_vars_cell                = {'MRep'};              % Inflation misreport
monpol_vars_label_cell          = {'Relative Underreport'};
monpol_vars_label_short_cell    = {'Mrep'};

factors1_cell_MRep              =   {'dBE'};     % Factor
factors1_label_cell_MRep        =   {'dBE'};     % Factor name

smpl_min_factors_vec =   [ones(1,1)*[2007 2]];         % Starting date for all of the factors
smpl_max_factors_vec =   [ones(1,1)*[2008 mtmp]];      % Finishing date for all of the factors.

VAR.irhor   = 20;                                      % Impulse Response Horizon
VAR.p       = 1;                                       % VAR lag length
VAR.SD_size = 1;                                       % Number of STD for the IMP RESP SHOCK

no_smpl_min_VAR = size(smpl_min_VAR_vec,1);
no_smpl_max_VAR = size(smpl_max_VAR_vec,1);
no_monpol_vars  = length(monpol_vars_cell);

ii_monpol  = 1;
ii_factors = 1;

ii{ii_monpol} = 1000*(ii_monpol-1)+1;
eval(['factors1_cell=factors1_cell_' monpol_vars_cell{ii_monpol} ';']);
eval(['factors1_label_cell=factors1_label_cell_' monpol_vars_cell{ii_monpol} ';']);

no_factors          =   length(factors1_cell);
no_factors1         =   length(factors1_cell);

for ii_smpl = [1]             
    switch ii_smpl
        case 1
            ii_smpl_min     =   1;
            ii_smpl_max     =   1;
        case 2
            ii_smpl_min     =   1;
            ii_smpl_max     =   2;
    end            

    smpl_min_VAR = smpl_min_VAR_vec(ii_smpl_min,:); % Get first date for sample in VAR

    % Factors sample starts minimum p periods after the VAR sample
    if smpl_min_factors_vec(ii_factors,1)>=smpl_min_VAR(1,1)+VAR.p/12
        smpl_min_FACTORS = smpl_min_factors_vec(ii_factors,:);
    else
        smpl_min_FACTORS(1,1) = smpl_min_VAR(1,1)+floor(VAR.p/12);
        smpl_min_FACTORS(1,2) = smpl_min_VAR(1,2)+(VAR.p-12*floor(VAR.p/12));
    end

    smpl_max_VAR = smpl_max_VAR_vec(ii_smpl_max,:);
    if (smpl_max_factors_vec(ii_factors,1)>smpl_max_VAR(1,1))
        smpl_max_FACTORS    =   smpl_max_VAR;
    else
        smpl_max_FACTORS = smpl_max_factors_vec(ii_factors,:);
    end

    % Print what is being calculated
    fprintf('\n\n#%3.0f\n',ii{ii_monpol});
    fprintf(['MONPOL: ' monpol_vars_cell{1,ii_monpol} ...
        '\nSAMPLE: ' num2str(smpl_min_VAR(1,1)) '-' num2str(smpl_max_VAR(1,1)) '\n']);                    

    % Find the dates in the sample
    VAR.smpl_min_VAR = find(and(DATASET_VAR.TSERIES(:,cell2mat(values(DATASET_VAR.MAP,{'YEAR'})))==smpl_min_VAR(1,1), ...
        DATASET_VAR.TSERIES(:,cell2mat(values(DATASET_VAR.MAP,{'MONTH'})))==smpl_min_VAR(1,2)));

    VAR.smpl_max_VAR = find(and(DATASET_VAR.TSERIES(:,cell2mat(values(DATASET_VAR.MAP,{'YEAR'})))==smpl_max_VAR(1,1), ...
        DATASET_VAR.TSERIES(:,cell2mat(values(DATASET_VAR.MAP,{'MONTH'})))==smpl_max_VAR(1,2)));

    VAR.smpl_max_VAR_factors = find(and(DATASET_VAR.TSERIES(:,cell2mat(values(DATASET_VAR.MAP,{'YEAR'})))==smpl_max_FACTORS(1,1), ...
        DATASET_VAR.TSERIES(:,cell2mat(values(DATASET_VAR.MAP,{'MONTH'})))==smpl_max_FACTORS(1,2)));        % Maximum place of factors in the VAR dataset
    
    VAR.smpl_min_FACTORS = find(and(DATASET_FACTORS.TSERIES(:,cell2mat(values(DATASET_FACTORS.MAP,{'YEAR'})))==smpl_min_FACTORS(1,1), ...
            DATASET_FACTORS.TSERIES(:,cell2mat(values(DATASET_FACTORS.MAP,{'MONTH'})))==smpl_min_FACTORS(1,2)));
    VAR.smpl_max_FACTORS = find(and(DATASET_FACTORS.TSERIES(:,cell2mat(values(DATASET_FACTORS.MAP,{'YEAR'})))==smpl_max_FACTORS(1,1), ...
            DATASET_FACTORS.TSERIES(:,cell2mat(values(DATASET_FACTORS.MAP,{'MONTH'})))==smpl_max_FACTORS(1,2)));

    %%% Select the variables in the VAR.
    % Start with policy indicator.
    VAR.select_vars             = {monpol_vars_cell{1,ii_monpol}};
    VAR.select_vars_label       = {monpol_vars_label_cell{1,ii_monpol}};
    VAR.select_vars_label_short = {monpol_vars_label_short_cell{1,ii_monpol}};
    ii_vars                     = 1;  % Size of Y vector in VAR.

    % Followed by spreads and IP cycle.
    VAR.select_vars             =  [VAR.select_vars,{'SP','IP'}];   % Add SP and IP. Order of vector is [MRep; SP; IP]. Will be reordered below given Cholesky.
    VAR.select_vars_label       =  [VAR.select_vars_label,{'Sovereign Spreads','Economic Activity'}];                 
    VAR.select_vars_label_short =  [VAR.select_vars_label_short,{'SP','IP'}];
    
    ii_vars        = ii_vars+2; % reorder in decreasing order of exogeneity (lower triangle)
    VAR.chol_order = [3 2 1];   % Cholesky ordering of the selected variables 
                                % [3=IP, 2=SP, 1=EcAct]
                                 
    VAR.select_factors          = {factors1_cell{1,ii_factors}};
    VAR.select_factors_label    = {factors1_label_cell{1,ii_factors}};

    % Run the VAR
    [VAR,VARChol,VARbs,VARCholbs] = BA_doVAR(VAR,DATASET_VAR,DATASET_FACTORS,nboot,clevel);

    nRow    =   3;
    nCol    =   2;
    Nfg     =   BB_plot_figure_sep(VAR, VARChol, VARbs, VARCholbs, nRow, nCol, ii_smpl);

    ii{ii_monpol} = ii{ii_monpol}+1;
end;
